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Abstract 

We investigate numerically spherically symmetric collapse of a scalar field 
in the semi-classical approximation. We first verify that our code reproduces 
the critical phenomena (the Choptuik effect) in the classical limit and black 
hole evaporation in the semi classical limit. We then investigate the effect of 
evaporation on the critical behavior. The introduction of the Planck length 
by the quantum theory suggests that the classical critical phenomena, which 
is based on a self similar structure, will disappear. Our results show that 
when quantum effects are not strong enough, critical behavior is observed. In 
the intermediate regime, evaporation is equivalent to a decrease of the initial 
amplitude. It doesn't change the echoing structure of near critical solutions. In 
the regime where black hole masses are low and the quantum effects are large, 
the semi classical approximation breaks down and has no physical meaning. 
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1 Introduction 



Spherical symmetric collapse of a scalar field was studied analytically by Christodoulou 
m who concluded that for near trivial initial data, the field disperses to future null in- 
finity leaving an empty space-time, and for "stronger" initial data, the field implosion 
forms a black-hole. This result was verified numerically by Goldwirth & Piran If 
we now characterize the "strength" of the initial data by some well defined parameter 
p which is monotonously growing with the "strength" (e.g. the initial amplitude of 
the field), we can expect that there will exist a critical value p* so that data with 
p < p* disperses while data with p > p* creates a black-hole. An intriguing question 
is what happens at p = p*. This question was investigated numerically by Choptuik 
0. His findings were 

• In the limit p —>■ p* the asymptotic behavior of the scalar field and the metric 
is universal - independent of the initial profile of the field. The field and metric 
also posses a symmetry - discrete self similarity (DSS) in the neighborhood of 
the critical solution's horizon. 

• The DSS is a periodical behavior at the origin with exponentially smaller periods 
in proper time: r i-^ re~^ where t = t* — t, t* being the time of formation of 
the critical black-hole, A ^ 3.44. 

• The mass of the black hole formed by marginally supercntical data has a scaling 
law M oc (p - p*y with 7 0.375. 

Perhaps one of the striking conclusions from these findings is that it is possible to 
produce a black-hole with an infinitesimal mass. The limiting case zero mass black- 
hole, the "choptuon" can be viewed as a collapsing ball of field energy for which the 
rate of collapse is exactly balanced by the energy loss by radiation, so that when the 
ball shrinks to zero radius all of its energy is radiated away 0]. The choptuon is 
actually a naked singularity, visible from null infinity although it is not generic 
and it is destroyed by an arbitrarily small perturbation. 

Hawking has shown that within the semi-classical approximation black holes 
evaporate thermally. The rate of evaporation depends on the black hole's mass via 
the relation 0: 
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with the proportionahty constant depending on the number of particle species that 
can be emitted by the black-hole at a given temperature. Integration of this relation 
yields a black-hole lifetime which is proportional to M^. This suggests that quantum 
effects will change the Choptuik effect since infinitesimally small black holes evaporate 
almost instantly. 

The semi-classical approach to gravitation stipulates that we can write the semi- 
classical Einstein equation 

G,u = 8vr(T0)^,, (2) 

where (T^) is the expectation value of the effective quantum energy-momentum ten- 
sor. This approach has been tested to give rise to evaporation of the black hole 
satisfying Eq. (|T]) in 0. 

We explore the evaporation of a black hole formed by collapse of a scalar field 
and the effect of the evaporation on the Choptuik phenomena. To this end, we 
investigate, numerically, the spherical collapse of a scalar field with the addition of an 
appropriate {T'^). We use {T'^) adapted from a 2D expectation value that reproduces 
the Hawking radiation ||^. This in turn forces us to solve the evolution equations. 
This is significantly harder in our coordinates then the usual method of solving the 
constraint equations used in and 0. In fact, this is the first time that a numerical 



solution using the evolution equations has been presented. Another problem with the 
above {T'^) is that the tensor is divergent at the origin. We renormalize the tensor 
to bypass this problem. This renormalization violates the Bianchi identity near the 
origin however the semi classical approximation is not valid there anyway. The code 
is tested with the classical problem to reproduce the Choptuik effect. After this we 
check black hole evaporation. Finally we examine the effect of evaporation on the 
Choptuik phenomena. 

The introduction of quantum effects via (T^) immediately introduces the Planck 
length scale denoted here by y/a. The addition of a length scale to the previously 
scale-less problem of the massless scalar field suggests that the critical phenomena 



present there will disappear W^ . However, our results show that the critical phenom- 
ena does not disappear within the regime of validity of the semi-classical approxima- 
tion. 

We discuss the problem of classical & semi-classical scalar field collapse in Sec. II 
& HI. Our numerical scheme is described in Sec. IV and our results in Sec. V. In Sec. 
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u=v 



r=0 




Figure 1: (a) The boundaries of the computational region. Initial values are given on the 
initial surface and boundary conditions on the origin. The shaded area is the integration 
area, (b) The grid defined by double null coordinates. The casual past of some point p is 
shaded 

VI we summarize our findings. Finally in Appendix A we derive several auxiliary 
quantities needed for the numerical scheme. 



2 Classical Scalar Field Collapse 

The characteristics of a massless scalar field are null so a double null coordinate 
system u, v is the "natural" system for this problem. In these coordinates the horizon 
(when it forms) is regular and there are no coordinate singularities. This enables us 
to study quantum effects that occur near the horizon. 
The metric is: 

ds^ = e^^dudv - r^dn\ (3) 

Where / and r are functions of u and v only. This definition is only unique up to a 
change of variables v' — fi{v),u' = f2{u). This ambiguity can be settled by a choice 
of the origin {r = 0) as u — v and by choosing either r or / on the initial surface. 
We choose r — on u — which is supposed to be sufficiently far in the past. This 
choice corresponds to an asymptotically fiat space-time. Regularity of the origin, 
with the above choice of coordinates, implies the following boundary conditions: 



r=0 
r=0 



—r 



r=0' 
r=0' 



(4a) 

(4b) 
(4c) 
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The Classical Einstein equations together with the field's equation of motion can 
be divided into two types: The dynamical equations 



V 



and the constraint equations: 



~^ 11.11 



-- ir u(l),v + r^v(l),u) 
r 



We have introduced here the auxiliary quantity r], 



V 



^2r „r,„ + ^e^^ 



=2/ 



m 



where m(M, v), is the mass inside a sphere of radius r{u, v) 



m 



^(l + 4r„r,e-2/). 



(5a) 
(5b) 

(5c) 

(6a) 
(6b) 

(7) 
(8) 



3 Black Hole Evaporation & The Semi-Classical 
Approximation 

We add to the classical equations an effective stress energy tensor {T'^) that describes 
the evaporation. Following we turn to 2D theories to obtain a tensor embodying 
Hawking radiation and adapt it to 4D. The expectation value of the 2D quantum 
energy- momentum tensor |T^, in double- null coordinates, is divided by 47rr^ to mimic 
the 4D radial dependence to give: 



a 



47rr^ 



f,UU f^u f,uv 

f,UV f,VV f^D 







(9) 



The constant a defines a length scale y/a which is of the order of the Planck scale. 
Beyond this scale the semi-classical approximation it is not valid. 

Unfortunately this tensor diverges at the origin. Note that for a null fiuid con- 
sidered by 0, the origin is always fiat and (T'^) vanishes identically there. Any 
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attempt to regularise this tensor must take into account that this tensor must be 
divergence-less (i.e. V(T'3) = 0) because of the Bianchi identity. A simple attempt 
to renormalize (T*^) by multiplying it by a general function Q{r) fails, we find that 
the Bianchi identity is violated unless Q = Const. Other attempts which were based 
on an introduction of angular terms failed as well and we had no choice but to ignore 
the Bianchi identity following H] and to multiply (T^) by 



1 + (a/r2)' 



The deviation from the Bianchi identity is significant only for rli^^Ja.^ for which the 
whole the semi classical approximation is no longer valid anyway. 
The resulting semi-classical dynamical equations are: 



(^47rQ0,„0,, - I) , (11a) 



= {\ -g) (-47r0,„0,^, + I) , (lib) 
4>,uv = -- {r,u4>,v + r,v4>,u) , (11c) 



and the constraint equations are: 



-r„„ + 2/,^r^ = 47rr02^ + gr (/,^^ - /2 ) , (12a) 
-r,uu + 2lur\u = ^nrcPl + Qr^l^u- fl). (12b) 



4 Numerical Scheme 

We choose m = as our initial surface and use the constraint equation Eq. (|12a| ) to 
determine the value of / given r and on this line. The integration proceeds, for 
each constant u, line from the origin v = u, where the boundary conditions Eq. (^) 
are imposed, up to fmax which is an arbitrarily chosen maximal value of v. In this 
way, the complete casual past of each point is known when the integration routine 
has to determine the field's value at that point. 

The straightforward numerical scheme becomes unstable once the gravitational 
effects become important (^ ~ 1). To stabilize the scheme, we have introduced the 
auxiliary quantity rj Eq. (|^ which is integrated on constant u lines. We use this value 
of 7] in the source for the evolution equations. It is worth noting that rj needs to be 
accurate only to first order since it is a source term. This modification stabilizes the 
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Figure 2: (a) The regular grid structure away from the origin. The source terms are 
calculated at Si and S2. rj is integrated along the dashed line and the derivatives are 
calculated at pi, p2 in the first iteration, (b) A typical grid arrangement near the origin. 
The point p is interpolated. Notice that du ^ dv (e.g. cell 1-2-2'-!') and that the first 
rectangular cell might differ from the rest of the cells (e.g O-l-l'-p as compared to 1-2-2'- 
1'). Points -1 and -1' are "fake" points calculated from the symmetry conditions. 

scheme. This is the first numerical scheme to solve this problem using the evolution 
equations in double null coordinates. 

A grid refinement algorithm was used in order to increase resolution at specific 
places such as the near-critical solutions and at the evaporation stage of the black- 
hole. 

We take m = to be the initial surface. On this surface, we specify an initial 
value for the scalar field 0. Then we choose r = |f on the initial surface and solve 
Eq. (|12a|) for /. The quantum terms should be negligible on the initial surface, so 
they were ignored while solving for /. We later verify that these terms are indeed 
negligible. The initial value equation for / is 

/,„ = 27T — (j)l = 2'KV(f)% (13) 

with /(O, 0) = as a boundary condition. Eq. (|T^) is integrated using a fourth order 
Runge-Kutta algorithm. 

Once the values of the fields on ti = are known the integration proceeds as 
follows: (all references to grid points and cells are illustrated in Fig. (^) 

The first step involves a triangular cell on the origin. Here we utilize the boundary 
conditions Eq. (^. Because space-time is fiat near the origin, the r direction is (1,-1) 
when taken as a vector in the {v,u) plane. Thus we can approximate ^ to order 
(Ar)^ in the current step by using the value of at the two last steps. 

• Using the value of 0i' and (ps" we obtain the value of 0o- 
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• Once the value of 0o is known, (p^^ and (p^uv are calculated on the origin between 
the points and 0' using also 0_i = 0p (again using Eq. and 0o')- 

• Now we calculate the value of f^uv = vl'^ = 47r/30^^ and using the known /o' 
and f-i = fp we obtain /q 

• Finally we calculate f^^ at the same point as and (p^^v and from them rj 
and using Eq. (^) and Eq. ( PT| ) (see appendix A for derivation of these 
equations) . 

All interpolations are done using a third degree polynomial such that the point be- 
ing interpolated always has two known points on each side (e.g. the values z_i', zqi, zy 
and Z21 are used to interpolate Zp). 

Away from the origin we have rectangular cells. For each step in the v direction 
three points are known (from previous steps) and the fourth (the future-most) is 
calculated (e.g ^2,^2' and z^i are known and z^ is calculated). The integration is a 
two step iterative procedure. The basic principle is that the second derivative operator 
can be discretesized as follows: 



Z^- Z2- Zy + Z2' 

AuAv 



+ 0{{AuAv) 



(14) 



The source terms of the equations are split into two - the gravitational term involving 
1] and the field term. The gravitational term is integrated along a constant u line 
which runs through the point at which the source is calculated (in between of the 
grid-lines). The field term is calculated directly. 

1^* Iteration The fields and their derivatives at point S2 are approximated using the 
known values at 2,2' and 3': 



S2 



S2 



S2 



1/ 

-{Z2 + Z2' + Zy) 



Av 



{Zy - Z2') 



Z2' 



7] is integrated using rj^ calculated for the previous point. Eq. (|T 

' SI 

first approximation for z^. 



(15) 
(16) 
(17) 

yields a 
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2"*^ Iteration A second order approximation for the fields and their derivatives can 
now be obtained using calculated in the first iteration: 

1 



4 



S2 



S2 



Z2 + ^3 + Z2' + 



{Zy - Z2' + ^3 - ^2), 



2/\v 

1 

2Au 



Z2 - Z2' + 23 - 2:3/ 



An improved approximation for r] is given by: 

1 



pi 



S2 



(18) 
(19) 
(20) 

(21) 



This iterative explicit scheme gives us a second order accuracy. 

To study the last stages of collapse, we utilize a cell refinement algorithm that 
increases the resolution of the integration in the areas of interest. Various multi grid 
or adaptive mesh |]13[ methods have been used to study critical collapse (cf. and 
0) but the increased complexity of these schemes was unnecessary for our purpose. 
We chose to half the value of Am at preset intervals, keeping the value of Av constant. 
This scheme has the advantage that we can reach an asymptotic value of u which is 
independent of the value of fmax- Using it we can study structure away from the 
origin. The price is that the first rectangular cell can be different from the others (see 
Fig. p|b). This scheme also sacrifices v resolution. Both problems are not critical. 



5 Results 

To test our code we first ran it on the classical (i/a = 0) problem. We find a 
critical initial amplitude p*. For initial data with amplitudes smaller then p* the 
field disperses leaving a fiat space-time. For initial amplitudes above p* the field 
collapses to form a black hole. For amplitudes near p* the solution exhibits self 
similar oscillations. A log-log plot of Mbh as a function of p — p* can be seen in Fig. 
d^). Although the critical exponent we found, 7 = 0.409 is larger then the usually 
quoted number 7 ^ 0.375 the power law dependence of Mbh on p—p* is evident. Thus 
our code reproduces the classical behavior. The discrepancy in the value of 7 reflects 
the accuracy of the code. We also present a contour plot of r showing the apparent 
horizon, the singularity and the event horizon (defined as the last ray to avoid the 
singularity). For the classical case, the apparent horizon is always inside the event 
horizon - any photon which starts falling into the black hole will never escape. 
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Figure 3: Classical results: (a) A log-log plot of Afbh vs. p — p* . The slope, 0.409, is 
higher then the usually quoted result of 0.375. Both the slope and p* were found through 
least squares fitting, (b) A contour plot of r values in the case of formation of a black hole. 
The event horizon, apparent horizon & the singularity are marked. The apparent horizon 
is always beyond (larger u) the event horizon. 



The introduction of (T*^) creates an effective outgoing flux corresponding to Hawk- 
ing radiation. This flux produces an evaporation satisfying Eq. (|I]) as verified in 
for black holes created by a null fluid. It is sometimes hard to distinguish the evap- 
oration from the energy momentum of the scalar field which is reflected from the 
origin. These reflections and the tail of the incoming scalar field can mask the quan- 
tum flux. To see a distinct black hole evaporation we need to create a situation where 
the only energy flux crossing the horizon is the outgoing Hawking radiation. This 
can be accomplished by looking at initial conditions that create a black hole before 
a significant fraction of the scalar field's energy is reflected from the origin. We also 
look for conditions producing an apparent horizon after most of the field has collapsed 
(large v) so that the tail of the field is also inside the horizon. 

Fig. (^ depicts the geometry resulting from the evolution of a scalar field with 
such initial conditions. In the r contour plot Fig. (^) The appearance of the apparent 
horizon before (lower u) the event horizon is apparent. Null trajectories which begin 
to curve back towards the singularity but then turn around and escape to infinity are 
shown in the r{v) plot Fig. (|^). The evaporation of the black hole is indicated by 
the decreasing radius of the outer apparent horizon which implies a negative energy 
flux through the apparent horizon. 

It is difficult to resolve the evaporation in the u, v coordinates which constitute 
the inertial frame of an observer at rest on the origin. In these coordinates, the 
evaporation is contained within a tiny u lapse. We therefore utilize the redshift of 
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Figure 4: (a) Contour of r in the evaporation region. The field is centered around v = 20, 
to the left of it the space time is flat (r oc w — u). The dotted hne describes r.t, = - the 
apparent horizon. The event horizon is at the top of the graph (u = 3.528). (b) trajectories 
of r{v) for different values of u. These null trajectories curve back towards the singularity 
but then turn around and escape to infinity. This is a clear indication for mass decrease. 
The decrease in radius of the apparent horizon is evident, (c) The black hole's mass as 
function of w. A 25% decrease in mass is evident. 
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Figure 5: Mhh as function of ^/a and p — p* Darker colors indicate lower mass. As ^/a 
increases Mbh decreases until finally no black hole is formed. The empty region is where 
no black holes formed. The final decrease in Mbh is too rapid to be resolved in this graph 

outgoing rays (as in HI) which imphes that diverges as the event horizon is 

approached. We half the u step Am in order to always have: 

r,„AM _ < Ar, (22) 

where Ar is the value of the Ihs. of Eq. (|2^) for the first u step. Using this scheme 
we managed to resolve a 25% decrease in the black hole mass due to evaporation. 

In our search for the scaling law which characterizes the Choptuik phenomena, 
we have to stop when the resulting black hole mass approaches ^/a. Nevertheless, we 
can still take a <^ 1 and try to get a scaling law down to the minimal meaningful 
mass. 

In Fig. (^ we show a plot of Mbh as a function of p—p* and ^/a. For small ^/a the 
evaporation is negligible and we recover the power law dependence of Mbh on p — p*. 
As we increase y/a the effects of mass loss is apparent. For large values of y/a black 
holes do not form where they would have classically. Thus, increasing a is effectively 
equivalent to decreasing p. This occurs because of mass loss during the evolution. 
These losses by evaporation can even change supercritical classical initial data into a 
subcritical solution and inhibit the formation of a black hole. The decrease in mass 
with increasing ^/a is too rapid to be resolved by the coarse grid of Fig. In Figs, 
d^) and (1^ we examine more closely the changes in the evolution due to increasing 
y/a. In Fig. (^ we show, Mbh as a function of y/a for three different values of p — p*. 
A power law dependence of Mbh on ^/a with an exponent which is dependent on the 
value of p — p* can be seen. The power law depicted in Fig. (|]) is broken before the 
semi-classical limit is reached by a rapid fall in Mbh- This feature is also present in 
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Figure 6: Mbh as function of a* — a for three different p. a* was determined through 
fitting. The values of a* are 3.8 x 10"*, 2.8 x IQ-^ and 2.7 x IQ-^^ for graphs (a), (b) and 
(c) respectively. 




Figure 7: f and on origin for fixed initial data (fixed p) and different ^/a. The different 
curves are offset vertically in order to enhance readability. Increasing ^/a increases the 
number of echoes and is thus equivalent to decreasing p 
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the classical {a = 0) case indicating that we have reached our limiting resolution. Fig. 
(J^) depicts and / on the origin for various values of y/a with a fixed value of p. We 
see that the effect of increasing ^/a is practically similar to that of decreasing p. The 
general echoing structure of near critical solutions is conserved and the decrease in 
mass affects only the last echo. This is in agreement with the power law dependence 
of Mbh on y/a since each echo is exponentially smaller then the preceding one in mass. 



6 Summary 

The numerical solution of spherical symmetric semi-classical scalar field collapse en- 
abled us to study the effects of black hole evaporation on the critical phenomena 
present in classical collapse (the Choptuik effect). The code was tested to reproduce 
the critical phenomena in the classical limit and black hole evaporation in the semi 
classical approximation. 

The addition of quantum effects and the subsequent introduction of the Planck 
scale ^/a changed the behavior of the solutions. For initial data found to form a black 
hole classically {a = 0) the introduction of ^/a caused a decrease in the resulting 
black hole mass Mbh • For large enough ^/a, no black hole formed. This behavior 
is expected due to mass loss through evaporation. A detailed investigation of the 
relations between the parameters revealed that for a fixed initial data (fixed p) there 
is a power law dependence of Mbh on y/a. The exponent of the power law depends on 



p. Unlike some earlier expectations ||TT| we find that the introduction of the Planck 



scale does not destroy altogether the DSS behavior observed classically. 



A Auxiliary Quantities 

In order to carry out the numerical integration of the equations ( pJ]) we need to 
calculate some auxiliary quantities, mainly get an equation for rj^^ and get the value 
of 7] and on the origin. 

Using Eq. (|]) we can get an equation for m from the constraint equation Eq. ( [L2a| ) 

/ 2m \ 

2r„m,, = An[l-—jr^(l)l (23) 
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using this we can calculate rj on the origin (which is always classic) 



m{u, v) 



271 1 



2m\ 



so by a change of variable from f to r at constant u we have 



m (m, r{u, v)) 



r ( 2m\ r'^ 
27r 1 I — 



J dr' 



now, to the lowest order in r, we have (since — oc o (r^)) 



Jo 



m = 27r^ / r''^dr' 



(24) 



(25) 



(26) 



where all the terms out of the integral are to be evaluated at r = 0. Using r\ = 
on the origin, we have, for small r 



m 



and so, on the origin, using Eq. (^ 



in 



to get an equation for r/ ^,, using, again, Eq. (|^ we have 



=2/ 



'n,v 



2m f y — 3m— + m^] 
r ' ) 



with Eq. d^) and the quantum corrections from Eq. ( p3|) the equation is 



2/,. - 3— h + — 2i]r 



(27) 



(28) 



(29) 



(30) 



+ Q 



=2/ 



— 2rir 



X 



f ,VV f V I 



We can find the value of the rhs. of Eq. ( pOD on the origin by expanding all the 
fields in powers of f — m. We then obtain the expansion coefficients (which depend 
on u) by substituting back into the Einstein equations Eq. (|n|). This is a rather 
tedious process but can be easily automated using symbolic math programs such as 
Mathematica & Maple. By comparison of coefficients we deduce that: 



(31) 



(all quantities are evaluated at the origin). 
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